pacman::p_load(sf, sfdep, tmap, tidyverse, DT)In-class Exercise 6
In-class
In this exercise we learn how to work with Spatial Autocorrelation
1. Setup
hunan <- st_read(dsn = "data/geospatial", layer = "Hunan")Reading layer `Hunan' from data source
`/Users/jeffery/Projects/Y4S1/IS415/S0methingSimple/IS415-GAA/in-class_ex/ex06/data/geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 88 features and 7 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 108.7831 ymin: 24.6342 xmax: 114.2544 ymax: 30.12812
Geodetic CRS: WGS 84
hunan2012 <- read_csv("data/aspatial/Hunan_2012.csv")hunan <- left_join(hunan,hunan2012) %>%
select(1:3,7,15,16,31,32)# nb: A neighbor list object created by st_neighbor
# style: "W" for standardized weights
# .before = 1: Insert into the front
# Returns a table with the neighbors that can be viewed
wm_q <- hunan %>%
mutate(nb = st_contiguity(geometry),
wt = st_weights(nb, style = "W"), .before = 1)
datatable(wm_q)moranI <- global_moran(wm_q$GDPPC, wm_q$nb, wm_q$wt)
glimpse(moranI)List of 2
$ I: num 0.301
$ K: num 7.64
K is average neighbours found
# Performs a basic test
global_moran_test(wm_q$GDPPC, wm_q$nb, wm_q$wt)
Moran I test under randomisation
data: x
weights: listw
Moran I statistic standard deviate = 4.7351, p-value = 1.095e-06
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.300749970 -0.011494253 0.004348351
p-value: < 0.05, enough statistical evidence we can reject null hypothesis, does not conform to random distribution (95% sure) - Fail to reject null hypothesis if greater than 0.05 no point proceeding Moran I: Greater than 0 suggest clustering in the data, but is a relatively low clustering
# Always best practice to set seed before simulation
set.seed(1234)Global Moran I Permutation
global_moran_perm(wm_q$GDPPC, wm_q$nb, wm_q$wt, nsim = 99)
Monte-Carlo simulation of Moran I
data: x
weights: listw
number of simulations + 1: 100
statistic = 0.30075, observed rank = 100, p-value < 2.2e-16
alternative hypothesis: two.sided
p-value: Is smaller Moran I: Is the same
Local Moran I Permutation
lisa <- wm_q %>%
mutate(local_moran = local_moran(GDPPC, nb, wt, nsim = 99), .before = 1) %>%
unnest(local_moran) # Turn it into a table form
datatable(lisa)ii: local moran i p_ii: p-value with base method p_ii_sim: Based on simulation mean: label hotspots median: Use if there is significant no. of positive or negative skew Plot out skewness to see if you want to use median
tmap_mode("plot")
pii_m <- tm_shape(lisa) +
tm_fill("p_ii") +
tm_borders(alpha=0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Local Moran's p-value", main.title.size = 1)
ii_m <- tm_shape(lisa) +
tm_fill("ii") +
tm_borders(alpha=0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Local Moran's I", main.title.size = 1)
tmap_arrange(pii_m, ii_m, ncol = 2)
lisa_sig <- lisa %>%
filter(p_ii < 0.05)
tm_shape(lisa) +
tm_polygons() +
tm_borders(alpha = 0.5) +
tm_shape(lisa_sig) +
tm_fill("median") +
tm_borders(alpha = 0.4)
Hot Spot Cold Spot Analysis
# MUST Use distance inverse distance, Further are smaller
wm_idw <- hunan %>%
mutate(nb = st_contiguity(geometry),
wts = st_inverse_distance(nb, geometry, scale = 1, alpha = 1),
.before = 1)HCSA <- wm_idw %>%
mutate(local_Gi = local_gstar_perm(GDPPC, nb, wts, nsim = 99), .before = 1) %>%
unnest(local_Gi) # Turn it into a table form
datatable(HCSA)gi_star value:
LISA use cluster and outlier G*: Hot and cold spot
h_pii_m <- tm_shape(HCSA) +
tm_fill("p_sim") +
tm_borders(alpha=0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Local G* p-value", main.title.size = 1)
h_ii_m <- tm_shape(HCSA) +
tm_fill("gi_star") +
tm_borders(alpha=0.5) +
tm_view(set.zoom.limits = c(6,8)) +
tm_layout(main.title = "Local G* I", main.title.size = 1)
tmap_arrange(h_pii_m, h_ii_m, ncol = 2)
HCSA_sig <- HCSA %>%
filter(p_sim < 0.05)
tm_shape(HCSA) +
tm_polygons() +
tm_borders(alpha = 0.5) +
tm_shape(HCSA_sig) +
tm_fill("gi_star") +
tm_borders(alpha = 0.4)